
%%% Weighting the Evidence: A Rank-Dependent Model of Outdoor Recreation, June 2024
%%% This file calibrates, for each of the three utility models considered, the total number of choice occasions so that the number of expected trips and landings match the actual data for New Jersey in 2021

%Importing the different data files. 
Tmain=readtable('ReadyDataFluke2022.xlsx'); 
cost=importdata('trip_costs_NE.xlsx'); 
cost=cost.data(cost.textdata(:,2)=="NJ");
C=unique(cost)/1000; 
y=[Tmain.age Tmain.income_medium Tmain.income_high Tmain.education_college Tmain.education_graduate Tmain.avidity Tmain.male];
Y=unique(y,'rows'); 
CatchPerTrip=importdata('CatchPerTrip2021.mat');
SFreshaffle=CatchPerTrip.SFreshaffle2021;
Sreshaffle=CatchPerTrip.Sreshaffle2021;
BSBreshaffle=CatchPerTrip.BSBreshaffle2021;
clearvars SFreshaffle21 Sreshaffle21 BSBreshaffle21 Tmain
parameters_rdu_tk=importdata('parameters_rdu_tk.txt');  
parameters_rdu_prelec=importdata('parameters_rdu_prelec.txt'); 
sigma_rdu_tk=importdata('sigma_rdu_tk.txt');  
sigma_rdu_prelec=importdata('sigma_rdu_prelec.txt');

SFbaglimit=3; 
SFsizelimit=18; 
Sbaglimit=50;
Ssizelimit=9;
BSBbaglimit=10;
BSBsizelimit=12.5;

NETrips=importdata('AvgeHarvestTrip.xlsx');
[r7,c7]=find(contains(NETrips.textdata,'NJ'));
NETrips.data=NETrips.data(r7,:);
[r8,c8]=find(NETrips.data==2021);
NETrips.data=NETrips.data(r8,:);
TRIPS=round(NETrips.data(1,2)/1000);
cutoff_sf=NETrips.data(1,3);
bsb_k=NETrips.data(1,4);
scup_c=NETrips.data(1,5);
sf_r=NETrips.data(1,6);
bsb_r=NETrips.data(1,7);
clearvars NETrips

T=100; 
M=25000; 

[AggregateKeepSFnumber_eu,AggregateReleasedSFnumber_eu,AggregateKeepSnumber_eu,AggregateKeepBSBnumber_eu,AggregateReleasedBSBnumber_eu,AggregateKeepSFnumber_mixed,AggregateKeepSFnumber_rdu_tk,AggregateKeepSFnumber_rdu_prelec,AggregateReleasedSFnumber_mixed,AggregateReleaseSFnumber_rdu_tk,AggregateReleaseSFnumber_rdu_prelec,AggregateKeepSFnumberLinear,AggregateReleasedSFnumberLinear,AggregateKeepSnumber_mixed]=deal(zeros(M,T));  
[AggregateKeepSnumber_rdu_tk,AggregateKeepSnumber_rdu_prelec,AggregateKeepSnumberLinear,AggregateKeepBSBnumber_rdu_tk,AggregateKeepBSBnumber_rdu_prelec,AggregateReleasedBSBnumber_rdu_tk,AggregateReleasedBSBnumber_rdu_prelec,AggregateKeepBSBnumber_mixed,AggregateKeepBSBnumberLinear,AggregateReleasedBSBnumberLinear,AggregateReleasedBSBnumber_mixed,TotalWelfare_mixed,TotalWelfare_eu,TotalWelfare_rdu_tk,TotalWelfare_rdu_prelec,TotalWelfareLinear,NumberOftrips_eu,NumberOftrips_mixed,NumberOftrips_rdu_tk,NumberOftrips_rdu_prelec,NumberOftripsLinear]=deal(zeros(M,T));
[vartheta_eu,vartheta_rdu_tk,vartheta_rdu_prelec,varthetaLinear]=deal(zeros(T,7)); 

tic;  
rng(0, 'combRecursive'); 

N=30; 

u_eu=zeros(T,N); 
[v_eu,opt_out_rdu_tk,opt_out_rdu_prelec,opt_out_eu,probTripRDU_tk,probTripEU,probTripRDU_prelec,v_rdu_tk,v_rdu_prelec,ReleasedSFnumber,ReleasedSnumber,ReleasedBSBnumber,BSBnumber,Snumber,SFnumber]=deal(zeros(M,T)); 
[SFkept,BSBkept,BSBreleased,SFkeptNumbers,SFreleased,Scatch]=deal(zeros(N,T)); 

myObj=@(p)pstar(M,N,SFreshaffle(:,:,1),p,cutoff_sf,SFbaglimit);
options=optimset('MaxFunEvals',50000,'TolX',1e-15,'MaxIter',10000, 'Display','on');
[p_cutoff_SF,fval,exitflag,output]=fzero(myObj,0.9,options);

par_rdu_tk=mvnrnd(parameters_rdu_tk,sigma_rdu_tk,T);
par_rdu_prelec=mvnrnd(parameters_rdu_prelec,sigma_rdu_prelec,1.2*T);
par_rdu_prelec(par_rdu_prelec(:,16)<0,:)=[]; 
par_eu=par_rdu_tk;
par_eu(:,length(parameters_rdu_tk))=1;

nbinR=zeros(T,1);
nbinP=zeros(T,1);
for t=1:T
     x=reshape(SFreshaffle(:,:,t),[],1);
     keep=zeros(1,length(x));
for j=1:length(x)
z=rand(1,x(j));
keep(j)=sum((z>=p_cutoff_SF),'all');
end
pd=fitdist(keep','nbin');
nbinR(t)=pd.R;
nbinP(t)=pd.P;
end

Cost_data=[datasample(C,M) datasample(C,M)];
Indcharac_data=datasample(Y,M,1);

%1. RANK-DEPENDENT UTILITY TVERSKY & KAHNEMAN
  k=0;
  while (k<M) && (mean(sum(NumberOftrips_rdu_tk,1))<=TRIPS)     
    
    k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end           
    end 
    
    vartheta_rdu_tk(t,:)=[par_rdu_tk(t,9),par_rdu_tk(t,10),par_rdu_tk(t,11),par_rdu_tk(t,12),par_rdu_tk(t,13),par_rdu_tk(t,14), par_rdu_tk(t,15)];
    v_rdu_tk(k,t)=RDU_tk(par_rdu_tk(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac);
    opt_out_rdu_tk(k,t)=par_rdu_tk(t,8)+IndvCharac*vartheta_rdu_tk(t,:)'+par_rdu_tk(t,7)*CostCharacOpt_out;
    probTripRDU_tk(k,t)=exp(v_rdu_tk(k,t))/(exp(v_rdu_tk(k,t))+exp(opt_out_rdu_tk(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*SFnumber(k,t);
    AggregateReleaseSFnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*sf_r;
    AggregateKeepBSBnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*bsb_k;
    AggregateReleasedBSBnumber_rdu_tk(k,t)=probTripRDU_tk(k,t)*bsb_r;  
    TotalWelfare_rdu_tk(k,t)=log(exp(v_rdu_tk(k,t))+exp(opt_out_rdu_tk(k,t)))/(-par_rdu_tk(t,7)/1000);
    NumberOftrips_rdu_tk(k,t)=probTripRDU_tk(k,t);
    end
  end

k1=k;
  
%2. RANK-DEPENDENT UTILITY PRELEC
  k=0;
  while (k<M) && (mean(sum(NumberOftrips_rdu_prelec,1))<=TRIPS)     
    
    k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N   
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
    end 
    
    vartheta_rdu_prelec(t,:)=[par_rdu_prelec(t,9),par_rdu_prelec(t,10),par_rdu_prelec(t,11),par_rdu_prelec(t,12),par_rdu_prelec(t,13),par_rdu_prelec(t,14), par_rdu_prelec(t,15)];
    
    v_rdu_prelec(k,t)=RDU_prelec(par_rdu_prelec(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac);
    opt_out_rdu_prelec(k,t)=par_rdu_prelec(t,8)+IndvCharac*vartheta_rdu_prelec(t,:)'+par_rdu_prelec(t,7)*CostCharacOpt_out;
    probTripRDU_prelec(k,t)=exp(v_rdu_prelec(k,t))/(exp(v_rdu_prelec(k,t))+exp(opt_out_rdu_prelec(k,t)));  
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*SFnumber(k,t);
    AggregateReleaseSFnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*sf_r;
    AggregateKeepBSBnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*bsb_k;
    AggregateReleasedBSBnumber_rdu_prelec(k,t)=probTripRDU_prelec(k,t)*bsb_r;  
    TotalWelfare_rdu_prelec(k,t)=log(exp(v_rdu_prelec(k,t))+exp(opt_out_rdu_prelec(k,t)))/(-par_rdu_prelec(t,7)/1000);
    NumberOftrips_rdu_prelec(k,t)=probTripRDU_prelec(k,t);
    end
  end

k2=k;

  
%3. EXPECTED UTILITY
    k=0;
  while (k<M) && (mean(sum(NumberOftrips_eu,1))<=TRIPS)
    k=k+1;
    IndvCharac=Indcharac_data(k,:); 
    CostCharac=Cost_data(k,1); 
    CostCharacOpt_out=Cost_data(k,2); 
    
    for t=1:T
    SFEncounters=SFreshaffle(k,:,t); 
    
    for j=1:N  
     if SFEncounters(j)>0
     z=rand(1,SFEncounters(j));
     SFlegal=(z>=p_cutoff_SF);
     SFkept(j,t)=min(SFbaglimit,sum(SFlegal,'all'));
     SFreleased(j,t)=SFEncounters(j)-SFkept(j,t);
     else
     SFkept(j,t)=0;
     SFreleased(j,t)=0;
     end
             
     u_eu(t,j)=(1-exp(-par_eu(t,6)*(par_eu(t,2)*SFkept(j,t)/100+par_eu(t,1)*sf_r/100+par_eu(t,4)*bsb_k/100+par_eu(t,3)*bsb_r/100++par_eu(t,5)*scup_c/100)))/par_eu(t,6);   
    end
   
    vartheta_eu(t,:)=[par_eu(t,9),par_eu(t,10),par_eu(t,11),par_eu(t,12),par_eu(t,13),par_eu(t,14), par_eu(t,15)];
    
    v_eu(k,t)=RDU_tk(par_eu(t,:),SFbaglimit,nbinR(t),nbinP(t),sf_r,bsb_k,bsb_r,scup_c,CostCharac);
    opt_out_eu(k,t)=par_eu(t,8)+IndvCharac*vartheta_eu(t,:)'+par_eu(t,7)*CostCharacOpt_out;
    probTripEU(k,t)=exp(v_eu(k,t))/(exp(v_eu(k,t))+exp(opt_out_eu(k,t)));
    SFnumber(k,t)=mean(SFkept(:,t));
    ReleasedSFnumber(k,t)=mean(SFreleased(:,t));
    
    AggregateKeepSFnumber_eu(k,t)=probTripEU(k,t)*SFnumber(k,t);
    AggregateReleasedSFnumber_eu(k,t)=probTripEU(k,t)*sf_r;
    AggregateKeepBSBnumber_eu(k,t)=probTripEU(k,t)*bsb_k;
    AggregateReleasedBSBnumber_eu(k,t)=probTripEU(k,t)*bsb_r;
    TotalWelfare_eu(k,t)=log(exp(v_eu(k,t))+exp(opt_out_eu(k,t)))/(-par_eu(t,7)/1000); 
    NumberOftrips_eu(k,t)=probTripEU(k,t);
    end
  end

k3=k;
  
varNames={'k_rdu_kt','k_rdu_prelec','k_EU'};
writetable(table(k1,k2,k3,'VariableNames',varNames),'k_NJ_2021.xls','Sheet',1);
save('p_cutoff_SF.txt','p_cutoff_SF','-ascii');

toc;
save('CalibratedSeasonNJ2021.mat');



